# Table D4: Covid-19 Unemployment Effects on Asian Hate Crimes: Nonlinear Models

rm(list = ls())

setwd('/path/to/replication/')

library(data.table)
library(MASS)
library(pscl)
library(qte)
library(xtable)

load('./data/panel_month_dummies.RData')

# Negative binomial model in DiD framework
# att
nb_b_c <- glm.nb(asian_hc ~ 1, data=panel[covid2==0 & losing_income_d==0])
nb_b_t <- glm.nb(asian_hc ~ 1, data=panel[covid2==0 & losing_income_d==1])
nb_a_c <- glm.nb(asian_hc ~ 1, data=panel[covid2==1 & losing_income_d==0])
nb_a_t <- glm.nb(asian_hc ~ 1, data=panel[covid2==1 & losing_income_d==1])

obs_atet <- mean((predict(nb_a_t, panel[covid2==1 & losing_income_d==1], type = 'response') -
                    predict(nb_b_t, panel[covid2==1 & losing_income_d==1], type = 'response')) -
                   (predict(nb_a_c, panel[covid2==1 & losing_income_d==1], type = 'response') -
                      predict(nb_b_c, panel[covid2==1 & losing_income_d==1], type = 'response')))

# p-value via randomization inference
if (!file.exists('./output/neg_binomial_ri.RData')) {
  source('./code/tableD4areg.R')
}

rr <- get(load('./output/neg_binomial_ri.RData'))
permutation_atet <- unlist(rr)
p_value <- (length(which(permutation_atet < -obs_atet)) + length(which(permutation_atet > obs_atet)))/length(permutation_atet)

nb <- data.frame(model = 'Negative binomial model in DiD framework',
           estimate = obs_atet,
           p_value = p_value,
           obs = nrow(panel)) 

# Zero-inflated negative binomial model in DiD framework
# att
panelu <- panel[!is.na(share_asian_immigrant)]

nbzi_b_c <- zeroinfl(asian_hc ~ 1 | share_asian_immigrant,
                           data=panelu[covid2==0 & losing_income_d==0], dist = 'negbin', link = 'logit')
nbzi_b_t <- zeroinfl(asian_hc ~ 1 | share_asian_immigrant,
                           data=panelu[covid2==0 & losing_income_d==1], dist = 'negbin', link = 'logit')
nbzi_a_c <- zeroinfl(asian_hc ~ 1 | share_asian_immigrant,
                           data=panelu[covid2==1 & losing_income_d==0], dist = 'negbin', link = 'logit')
nbzi_a_t <- zeroinfl(asian_hc ~ 1 | share_asian_immigrant,
                           data=panelu[covid2==1 & losing_income_d==1], dist = 'negbin', link = 'logit')

atet <- mean((predict(nbzi_a_t, panelu[covid2==1 & losing_income_d==1], type = 'response') -
                predict(nbzi_b_t, panelu[covid2==1 & losing_income_d==1], type = 'response')) -
               (predict(nbzi_a_c, panelu[covid2==1 & losing_income_d==1], type = 'response') -
                  predict(nbzi_b_c, panelu[covid2==1 & losing_income_d==1], type = 'response')))
obs_atet <- atet

# p-value via randomization inference
if (!file.exists('./output/zi_neg_binomial_ri_tmp.RData')) {
  source('./code/tableD4breg.R')
}

rr <- get(load('./output/zi_neg_binomial_ri_tmp.RData'))
permutation_atet <- unlist(rr)
p_value <- (length(which(permutation_atet < -obs_atet)) + length(which(permutation_atet > obs_atet)))/length(permutation_atet)

nbzi <- data.frame(model = 'Zero-inflated negative binomial model in DiD framework',
                 estimate = obs_atet,
                 p_value = p_value,
                 obs = nrow(panelu))

# Changes-in-Changes
panel_cic <- panel[,.(sum(asian_hc),max(losing_income_d)), by=.(covid2, panelvar)]
setnames(panel_cic, c('V1', 'V2'), c('asian_crimes', 'unemployment'))

# att
econ_cic <- CiC(asian_crimes ~ unemployment,
                           t=1, tmin=0, tname='covid2',
                           data=panel_cic, se = FALSE,
                           idname = 'panelvar'
)
obs_atet_cic <- econ_cic$ate

# p-value via randomization inference
# ri function
permute_qte <- function(permutation_data, n_reps) {
  results <- vector(mode='list', length=n_reps)
  for (r in seq(n_reps)) {
    print(paste("Rep", r))
    true_treatments <- data.table(permutation_data)[, .(is_treated=unique(unemployment)==1), by=panelvar]
    
    permutation_treatments <- data.table(
      is_treated_permutation = sample(true_treatments$is_treated),
      permutation_id = unique(permutation_data$panelvar)
    )
    
    permutation_data$is_treated_permutation <- NULL
    permutation_data <- merge(permutation_data, permutation_treatments, by.x='panelvar', by.y='permutation_id')
    
    econ_cic_count <- qte::CiC(asian_crimes ~ is_treated_permutation,
                               t=1, tmin=0, tname='covid2',
                               data=permutation_data, se = FALSE,
                               idname = 'panelvar'
    )
    
    results[[r]] <- econ_cic_count$ate
  }
  return(results)
}

set.seed(879)
rr_cic <- permute_qte(panel_cic, 500)

permutation_atet <- unlist(rr_cic)
p_value_cic <- (length(which(permutation_atet < -obs_atet_cic)) + length(which(permutation_atet > obs_atet_cic)))/length(permutation_atet)

cic <- data.frame(model = 'Changes-in-Changes framework',
                   estimate = obs_atet_cic,
                   p_value = p_value_cic,
                   obs = nrow(panel))


### table in article
print(xtable(rbind(nb, nbzi, cic),
       type = "latex",
       digits=4,
       caption = "Covid-19 Unemployment Effects on Asian Hate Crimes: Nonlinear Models"),
      caption.placement = "top", file = './tableD4.tex')

